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The discrete-dipole approximation (DDA) is a powerful method for calcu- 
lating absorption and scattering by targets that have sizes smaller than or 
comparable to the wavelength of the incident radiation. We present a new 
prescription - the Surface-Corrected Lattice Dispersion Relation (SCLDR) - 
for assigning the dipole polarizabilities that takes into account both target 
geometry and finite wavelength. We test the SCLDR in DDA calculations 
using spherical and ellipsoidal targets and show that for a fixed number 
of dipoles, the SCLDR prescription results in increased accuracy in the 
calculated cross sections for absorption and scattering. We discuss extension 
of the SCLDR prescription to irregular targets. © 2008 Optical Society of 
America 
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1. Introduction 

The discrete-dipole approximation (DDA) is a numerical technique for calculating 
scattering and absorption of electromagnetic radiation by targets with sizes smaller 
than or comparable to the incident wavelength. The method consists of approximating 
the target by an array of polarizable points (dipoles), assigning polarizabilities at these 
locations based on the physical properties of the target, and solving self-consistently 
for the polarization at each location in the presence of an incident radiation field. 
This procedure can yield arbitrarily accurate results as the number of dipoles used 
to approximate the target is increased. However, computational considerations limit 
the number of dipoles that can be used. Hence, methods for increasing the accuracy 
for a fixed number of dipoles are desirable. 

A key factor in determining the level of accuracy that can be reached for a 
given number of dipoles is the prescription for assigning dipole polarizabilities. In 
this work, we present a new polarizability prescription that takes into account both 
target geometry and the finite wavelength of incident radiation. We test this technique 
in calculations of absorption and scattering by spherical and ellipsoidal targets and 
show that for a fixed number of dipoles, it generally provides increased accuracy 
over previous methods. In Section El we discuss previous polarizability prescriptions 
and develop the new method. In Section |3] we present calculations testing the new 
prescription, and in Section |U we discuss our results. 
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2. Polarizability Prescriptions 

A fundamental requirement of the DDA is that the inter-dipole separation d be small 
compared to the wavelength of incident radiation, kd < 1, where k = u/c is the 
wavenumber in vacuo. Here we will assume the dipoles to be located on a cubic 
lattice with lattice constant d, as this facilitates use of fast-Fourier transform (FFT) 
techniques.^ 

The first implementations of the DDA^ used the so-called Clausius-Mossotti re- 
lation (CMR) to determine the dipole polarizabilities. In this procedure, the polariz- 
ability a is given as a function of the (complex) refractive index m as 

3d 3 fm 2 -l\ 

This approach is valid in the infinite wavelength limit of the DDA, kd — > 0. 

Draine 3 showed that for finite wavelengths, the optical theorem requires that the 
polarizabilities include a "radiative-reaction" correction of the form 

a ( nr ) 

a = 1- (2/3)i(a^ r )/d 3 )(kd) 3 ' ^ 

where is the "non-radiative" polarizability, that is, before any radiative-reaction 
correction is applied. Draine 3 used ocmr as the non-radiative polarizability. 

Based on analysis of an integral formulation of the scattering problem, Goedecke 
& O'BrierPand Hage & Greenberg^ suggested further corrections to the CMR polariz- 
ability of order (kd) 2 . Draine & GoodmarP studied electromagnetic wave propagation 
on an infinite lattice; they required that the lattice reproduce the dispersion relation 
of a continuum medium. In this "Lattice Dispersion Relation" (LDR) approach, the 
radiative-reaction correction emerges naturally, and the polarizability is given [to or- 
der (kd) 3 ] by 

" LDR = 1 + (a(°)/d 3 )[(&i + m 2 b 2 + m 2 b 3 S)(kd) 2 - (2/3)i(kd) 3 } ' (3) 

where = Qqmr is the polarizability in the limit kd — ► 0, b\ — —1.8915316, 
62 = 0.1648469 and 63 = —1.7700004, and S is a function of the propagation direction 
and polarization of the incident wave. S is given as 

S = E(%^) 2 , (4) 

3 

where a and e are the unit propagation and polarization vectors, respectively. Note 
that eq. (jlj) gives S = for waves propagating along any of the lattice axes. This 
method correctly accounts to 0[(kd) 3 } for the finite wavelength of incident radiation, 
and by construction, it accurately reproduces wave propagation in an infinite medium. 
Its primary limitation is that the accuracy in computing absorption cross-sections of 
finite targets (for a given number of dipoles) degrades rapidly as the imaginary part 
of the refractive index m becomes large (e.g., for Im(m) > 2). 
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A. Geometric Corrections: the Static Case 

Recently Rahmani, Chaumet & Bryant (RCB) proposed a new method for assigning 
the polarizabilities that takes into account the effects of target geometry on the local 
electric field at each dipole site. Consider a continuum target in a static, uniform 
applied field E°. At each location % in the target, the macroscopic electric field E™ is 
linearly related to E°: 

E™ = C^E (5) 

where C" 1 is a 3 x 3 tensor that will depend on location i, the global geometry of the 
target, and its (possibly nonuniform) composition. If we now represent the target by 
a dipole array, and require that the electric dipole moment Pj of dipole i be equal to 
d 3 times the macroscopic polarization density at location i, we obtain 

If on is the polarizability tensor of dipole i, then 

Pi = 

where —AijPj is the contribution to the electric field at location i due to dipole Pj 
at location j (this defines the 3x3 tensors Ay). Substituting (0) into (J7J) we obtain 

where the 3x3 tensors 

Aj = Ci-J2 A i3 (PC^d (9) 

can be evaluated (and easily inverted) if the Cj are known. 

The RCB approach requires that the tensors C first be obtained. For certain 
simple geometries, the Cj can be obtained analytically. For example, for homogeneous 
ellipsoids, infinite slabs, or infinite cylinders, the tensors C, can be expressed in the 
form 

C = 1 + (^) L ( 10 ) 

where L is a "depolarization tensor". For example, L = 1/3 for a homogeneous sphere. 

In the present work, we combine the LDR and RCB approaches in order to obtain 
a polarizability prescription that accounts both for finite wavelength and for local field 
corrections arising from target geometry. We adopt aRCB as the polarizability a<°> in 
the limit kd —>■ 0, and apply corrections up to 0[(kd) 3 ] based on the LDR. A further 
analysis of the electromagnetic dispersion relation of a non-cubic lattice^ called into 
question the value of the constant b 3 in eq. used by Draine & Goodman,© and 
found it instead to be undetermined by available constraints. Thus we include an 
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adjustable factor / whose value is chosen to optimize the behavior of the new method 
as discussed in the next section. The "Surface-Corrected Lattice Dispersion Relation" 
(SCLDR) polarizability is given by 

"SCLDR = ClRCB {1 + (a RCB /ci 3 )[(6i + m% + m %fS)(kd) 2 - (2/3)z(^) 3 ]}-\ (11) 

where 

/ = exp[-0.5Im(m) 2 ]. (12) 

In the next section, we test this new prescription in calculations of absorption and 
scattering by spherical and ellipsoidal targets. 

3. Sphere and Ellipsoid Calculations 

For a continuum target of volume V, the effective radius a cS = (317471-) 1 / 3 , the radius 
of a sphere of equal volume. The target is approximated by an array of N dipoles 
located on a cubic lattice, with the dipole locations selected by some criterion designed 
to approximate the shape of the original target. The inter-dipole spacing is then set 
to d = {V/Nf/\ 

For a given orientation of the dipole array relative to the incident wave, we 
calculate the cross sections C sca and C a b s for scattering and absorption, and the di- 
mensionless efficiency factors Q sca = C sc Jnal s , Q &hs = C ab s/Va 2 ff . 

To test the performance of the SCLDR polarizability prescription against pre- 
vious results, we performed a series of calculations using the DDA code DDSCAT,^ 
modified to permit use of the SCLDR polarizabilities. We computed Q sca and Q a b s 
for spherical targets with a range of refractive indexes and for a range of scattering 
parameters x = 27ia e s/X = ka c s, using three different approaches for assigning the 
dipole polarizabilities: LDR, RCB and SCLDR. Spherical targets were employed be- 
cause the exact optical properties can be readily calculated using Mie theory. We also 
performed a similar but more limited set of calculations for ellipsoidal targets. 

We tested the LDR, RCB and SCLDR prescriptions for a number of different 
refractive indexes in the region of the complex plane with Re(m) < 5 and Im(m) < 4. 
We determined that for refractive indexes with these ranges of real and imaginary 
parts, it was desirable for the SCLDR correction factor / to tend toward unity for 
Im(m) < 1 and to tend toward zero for Im(m) > 2. We chose the functional form of 
eq. (J12|) in order to reproduce this asymptotic behavior. 

Figures E and El show the results of calculations for spheres with refractive indices 
m = 1.33 + O.Oli and m = 5 + 4i, each approximated by an array of N = 7664 dipoles. 
Because the dipole array is not rotationally symmetric, Q sca and Q a b s calculated with 
the DDA depend in general on the target orientation; we perform calculations for 
12 orientations, and we show the average and range of the results. We calculate the 
fractional errors in Q sca and Q a bs by comparison with exact results obtained using 
Mie theory: 

& — Q{DD o<uT e) < 13 > 

Q(Mie) 

In previous workP^it was recommended that the DDA be used only when \m\kd < 
1, or a more stringent condition \m\kd < 0.5 if the DDA is to be used to calculate 
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the differential scattering cross section. In the present work we find that when the 
SCLDR polarizabilities are used, the fractional errors in Q sca and Q a bs are relatively 
insensitive to x provided \m\kd < 0.8, which we adopt as an operational validity 
criterion. Figures Q and El show results for values of x satisfying \m\kd < 0.8. 

From Figure^ it is clear that the LDR and SCLDR prescriptions provide approx- 
imately equal levels of accuracy in the \m\ ~ 1 regime, while the RCB prescription 
does not perform as well. Figure El shows that at the other extreme of Re(m) ^> 1 
and Im(m) ^> 1, the LDR approach results in large errors, especially in the calcu- 
lated absorption cross sections, while the RCB and SCLDR prescriptions perform 
approximately equally well. 

In Figures El and EJ we show the convergence behavior of the different polarizabil- 
ity prescriptions as the number of dipoles iV is increased for spherical targets with 
selected refractive indices; the refractive indices have been chosen to sample the re- 
gion of the complex plane discussed in the previous paragraphs. The SCLDR method 
performs comparably to or better than the RCB and LDR prescriptions throughout 
this region of the complex refractive index plane. This illustrates the advantage of 
the SCLDR approach over these previous techniques: it performs well not just for a 
small range of refractive indexes, but for the entire range we have sampled. 

Figures El and El extend the result shown in Figures El and |U to targets of a 
more general shape, specifically ellipsoids with approximately 1:2:3 axial ratios. For 
these targets, we have estimated the true values of Q scau and Qabs by assuming these 
to be linear functions of iV _1//3 , extrapolating to iV -1 / 3 — > for each polarizability 
prescription, and taking the average of the results from the different prescriptions. The 
close similarity of the results of these calculations to those shown in Figures El and El 
demonstrates that the SCLDR prescription provides the same benefits in calculations 
for ellipsoidal targets as for spheres, although we note that for ellipsoids with values 
of m with large imaginary parts [typically Im(m) > 1], the RCB prescription can 
provide improved accuracy in calculations of Q sca . 

For an isotropic material with refractive index m, the Clausius-Mossotti polariz- 
ability ttcMR has triply-degenerate eigenvalues «cmr = ( m2 — l)d 3 /4:7i. For the case of 
a 1:2:3 ellipsoid with refractive index m = 5 + 4i, we have calculated the eigenvalues 
a of ccscldr for \m\kd — > (for which case <*scldr — > «rcb) at each occupied lattice 
site. Figure [7| (left panel) shows the distribution of the fractional difference of the 
eigenvalues a from ocmr- The deviations tend to be appreciable (fractional difference 
exceeding ~20%) only near the surface. The left panel shows that the deviations ex- 
ceed 20% for 47% of the lattice sites for N = 90, but only 9% of the lattice sites when 
N = 43416. For this example the fraction of the eigenvalues deviating by >20% is 
~ 0.30(A r /1000)^ 1 / 3 for N > 500, approximately equal to the fraction of the dipoles 
located within a surface layer of thickness ~ 0.6d. 

The right panel in Figure [7| shows the eigenvalue deviations as a function of 
distance from the surface of the ellipsoid: the eigenvalues deviating from a CMR by 
more than ~ 20% are, as expected, exclusively associated with dipoles located within 
a distance d of the surface. 
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4. Conclusion 



We introduce a new DDA polarizability prescription - the Surface- Corrected Lattice 
Dispersion Relation (SCLDR). This technique builds on previous work, principally 
by Draine & GoodmarPand Rahmani, Chaumet & Bryant,^ to account properly for 
both finite wavelength and target geometry. We have tested the new polarizability 
prescription in calculations of absorption and scattering by spherical and ellipsoidal 
targets. These tests show that the SCLDR performs generally better than previous 
prescriptions which took account either of finite wavelength or of target geometry but 
not both. The SCLDR technique is most easily applicable to target shapes for which 
there exists an analytical solution to the electrostatic applied field problem, but it 
can be applied to any dielectric target (homogeneous or inhomogeneous, isotropic or 
anisotropic) provided that the electrostatic problem can at least be solved numerically 
to obtain the tensors Q (see eq. |3J). In such cases, it generally provides a significant 
increase in accuracy over previous methods, especially for highly absorptive materials. 
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Fig. 1. Comparison of scattering and absorption efficiency factors Q sca , and 
Q a bs computed for a pseudo-sphere of N = 7664 dipoles and refractive index 
m = 1.33 + O.Oli, averaged over 12 orientations, and using three different po- 
larizability prescriptions: Lattice Dispersion Relation (LDR); Rahmani et alP 
(RCB); and Surface- Corrected Lattice Dispersion Relation (SCLDR). The hor- 
izontal axis shows (top) \m\kd (the phase shift in radians within one lattice 
spacing) and (bottom) the scattering parameter x = ka. Error bars indicate 
the ranges of Q values obtained for the individual orientations. The top panel 
shows the results of Mie theory calculations; the lower panels show the frac- 
tional error in Q sca and Q a bs, respectively. The SCLDR and LDR prescriptions 
are clearly preferred over RCB for this case. 




Fig. 2. Same as Fig. ^ but for refractive index m = 5 + 4z. The SCLDR 
and RCB prescriptions are clearly preferred over the LDR for this case, with 
SCLDR being somewhat superior to RCB. 
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Fig. 3. Fractional error in Q a bs averaged over 12 orientations for spheres with 
different refractive indexes, as a function of iV -1 / 3 , where N is the number of 
dipoles, in the range 624-59278. Calculations are shown for the LDR, RCB, 
and SCLDR polarizability prescriptions; the symbolic scheme is the same as in 
Fig. Refractive indexes m, scattering parameters x = ka, and exact values 
of Qabs computed from Mie theory are shown in the left portion of each panel. 
The scattering parameters are chosen so that \m\kd ~ 0.8 (the approximate 
limit of applicability of the DDA) for the smallest number (N = 624) of 
dipoles. The convergence with increasing iV is quite smooth in all regions 
of the complex m-plane, with the exception of m = 3 + O.Oli. In almost 
every case shown, fractional errors < 2% (and often significantly lower) can be 
achieved for iV w 6000 dipoles. We find that for calculating Q a b s , the SCLDR 
is comparable or superior in accuracy to the LDR and RCB prescriptions 
throughout the region of m-space shown. 
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Fig. 4. Same as Fig. El except that fractional errors in Q sca are plotted. Again, 
the SCLDR prescription is comparable or superior to the LDR prescription for 
all values of m shown. While the SCLDR prescription is still comparable or 
superior to the RCB prescription for m values with small imaginary parts, the 
RCB prescription provides better accuracy in calculating Q SC3i for m values 
with large imaginary parts. 
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Fig. 5. Same as Fig. El but for ellipsoids with approximately 1:2:3 axial ratios. 
Fractional errors have been estimated based on comparison with an extrapo- 
lation of the convergence behavior of the three polarizability prescriptions, as 
described in Sectional Again the SCLDR prescription appears comparable or 
superior to the LDR and RCB prescriptions for calculating Q a b s throughout 
the region of the complex m-plane sampled, though the RCB prescription is 
slightly preferred for the cases of m = 3 + 4i and m = 5 + Ai. 
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Fig. 6. Same as Fig. El except that fractional errors in Q sca are plotted. As in 
Fig. El SCLDR is comparable to or superior to LDR for all values of m, and 
to RCB for values of m with small imaginary parts, while RCB is somewhat 
superior to SCLDR for values of m with large imaginary parts. 
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Fig. 7. Comparison of RCB and CMR polarizabilities. The left panel shows 
the distribution of polarizability eigenvalues for discrete dipole approximations 
to a 1:2:3 ellipsoid with m = 5+4i using N = 90, 688, 5456, and 43416 dipoles. 
The shaded region corresponds to a fractional difference of 20% or less; the 
fraction of the eigenvalues within this region varies from 53% for N = 90 
(3dxQdx9d axes) to 91% for N = 43416 (24dx48dx72d axes). The right panel 
shows the fractional difference between RCB and CMR polarizabilities versus 
the distance (in units of the shortest axis) from the ideal ellipsoidal surface 
used to define the target (all dipole locations are interior to this surface). As 
expected, the RCB polarizability reduces to the CMR polarizability for dipoles 
lying more than ~2d from the surface. 
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